Apêndice 2: Descrição Estatística
Mori, Danilo
31/10/2019
library(parallel)
library(doMC)
library(DHARMa)
library(broom.mixed)
library(gratia)
library(MuMIn)
library(mgcv)
library(bbmle)
library(lme4)
library(merTools)
library(magrittr)
library(gridExtra)
library(ggplot2); theme_set(theme_classic())
library(tidymv)
library(purrr)
library(plyr)
library(dplyr)1 - Variables
table 1: Variáveis utilizadas na descrição estatística
| code | description | class | range |
|---|---|---|---|
| n_nRef | number of SADs not refuted by the KS bootstrap test with critical p-value 0.05 | continuous | (0 ; 100) |
| diffS_mean | diffS = (S_MN - S_obs) / S_obs, diffS_mean = mean(diffS); S_MN = species richness estimated by MN | continuous | (-0.977 ; 4.78) |
| U_med | average of 10 replicates of the speciation rate estimated by MNEE | continuous | (8.860e-5 ; 4.221e-2) |
| p | proportion of tree cover | continuous | (0.013 ; 0.961) |
| MN | Neutral Model (EE = spatial explicit; EI = spatial implict) | category | 2 levels |
| d | mean dispersal distance (meters) | continuous | (0.456 ; 19.183) |
| k | proportion of propagules until the nearest neighborhood seq(0.99 : 0.05) | category | 20 levels |
| d_Lplot | mean dispersal distance / side of the sample area | continuous | (0.02 ; 0.192) |
| S_obs | observed species richness | integer | (26 ; 458) |
| Ntotal | number of individuals in the sample area | integer | (649 ; 12105) |
| SiteCode | sample area code | category | 103 levels |
Figure 1 Possible Predictor Variables
2 - Number of unrefuted SADs
Figura 2.1 número de SADs não refutadas ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.
Figura 2.2 número de SADs não refutadas ~ (d * MN|Site ~ p). Site está ordenado pelo valor de p.
2.1 GLMM binomial
Modelo cheio 1:
linear term
- ~ (p + p^2) * (d + d^2) * MN; if : var_dispersao = contiguous
- ~ ( p + p^2) * k * MN; if : var_dispersao = category
random term
- 1|SiteCode
- MN|SiteCode
- (var_dispersao + var_dispersao^2)*MN|SiteCode, if : var = dispersao : contiguous
# dados
df_md <- df_resultados %>%
select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>%
gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z)
# funcao ajuste dos modelos
f_md <- function(dados){
var_dispersao <- unique(dados$var_dispersao)
if(var_dispersao == "k"){
l_md <- vector("list",2)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"))
dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}else{
dados$value_dispersao <- as.numeric(dados$value_dispersao)
l_md <- vector("list",3)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"),
paste0(var_dispersao," (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}
return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio) <- NULL
l_nRef__modeloCheio <- do.call(c,l_nRef__modeloCheio)
#
f_warning <- function(glmm_object){
# glmm_object <- l_nRef__modeloCheio[[4]]
v_message <- glmm_object@optinfo$conv$lme4$messages %>%
as.character()
if(length(v_message)==0){v_message <- "OK"}
if(length(v_message)>1){v_message <- v_message[1] }
return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio,f_warning,.id="glmer") %>%
rename(warning_message = V1)Table 2.1 Modelos Cheios estimados e avisos de convergência
| glmer | warning_message |
|---|---|
| d_Lplot.z 1|Site | Model failed to converge with max|grad| = 0.00343548 (tol = 0.002, component 1) |
| d_Lplot.z MN|Site | Model failed to converge with max|grad| = 0.00836802 (tol = 0.002, component 1) |
| d_Lplot.z (d_Lplot.z+d_Lplot.z^2)*MN|Site | OK |
| d.z 1|Site | Model failed to converge with max|grad| = 40.1827 (tol = 0.002, component 1) |
| d.z MN|Site | Model failed to converge with max|grad| = 0.00746139 (tol = 0.002, component 1) |
| d.z (d.z+d.z^2)*MN|Site | OK |
| k 1|Site | OK |
| k MN|Site | OK |
allFit
i <- 1
names(l_nRef__modeloCheio[v_glmerUpdate])[i]
#
## update1:
md_allFit <- allFit(l_nRef__modeloCheio[v_glmerUpdate][[i]],maxfun = 1e5, parallel = 'multicore', ncpus = 3)1) d_Lplot 1|Site
7 optimizer(s) failed
2) d_Lplot MN|Site
7 optimizer(s) failed
3) d 1|Site
7 optimizer(s) failed
4) d MN|Site
7 optimizer(s) failed
Modelo cheio 2:
# dados
df_md <- df_resultados %>%
select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>%
gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z)
# funcao ajuste dos modelos
f_md <- function(dados){
var_dispersao <- unique(dados$var_dispersao)
if(var_dispersao == "k"){
l_md <- vector("list",2)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"))
dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}else{
dados$value_dispersao <- as.numeric(dados$value_dispersao)
l_md <- vector("list",4)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"),
paste0(var_dispersao," ",var_dispersao,"*MN|Site"),
paste0(var_dispersao,"^2 (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(value_dispersao * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[4]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}
return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio2 <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio2) <- NULL
l_nRef__modeloCheio2 <- do.call(c,l_nRef__modeloCheio2)
#
f_warning <- function(glmm_object){
v_message <- glmm_object@optinfo$conv$lme4$messages %>%
as.character()
if(length(v_message)==0){v_message <- "OK"}
if(length(v_message)>1){v_message <- v_message[1] }
return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio2,f_warning,.id="glmer") %>%
rename(warning_message = V1)Table 2.2 Modelos Cheios estimados e avisos de convergência
| glmer | warning_message |
|---|---|
| d_Lplot.z 1|Site | Model failed to converge with max|grad| = 0.00284637 (tol = 0.002, component 1) |
| d_Lplot.z MN|Site | OK |
| d_Lplot.z d_Lplot.z*MN|Site | OK |
| d_Lplot.z^2 (d_Lplot.z+d_Lplot.z^2)*MN|Site | OK |
| d.z 1|Site | Model failed to converge with max|grad| = 0.00369235 (tol = 0.002, component 1) |
| d.z MN|Site | OK |
| d.z d.z*MN|Site | OK |
| d.z^2 (d.z+d.z^2)*MN|Site | OK |
| k 1|Site | OK |
| k MN|Site | OK |
allFit
d_Lplot.z 1|Site
7 optimizer(s) failed
d.z 1|Site
7 optimizer(s) failed
Comparação de Modelos Cheios:
Tabela 2.3 Comparação baseada em AICc dos modelos cheios
| GLMM | dAICc | df | weight |
|---|---|---|---|
| d^2 (d+d^2)*MN|Site | 0.0 | 39 | 1 |
| d_Lplot^2 (d_Lplot+d_Lplot^2)*MN|Site | 274.6 | 39 | <0.001 |
| d d*MN|Site | 38427.3 | 22 | <0.001 |
| d_Lplot d_Lplot*MN|Site | 38586.5 | 22 | <0.001 |
| k MN|Site | 74141.2 | 123 | <0.001 |
| d MN|Site | 95710.0 | 15 | <0.001 |
| d_Lplot MN|Site | 96279.6 | 15 | <0.001 |
| k 1|Site | 106158.0 | 121 | <0.001 |
Tabela 2.4 Coeficiente de Determinação Condicional e Marginal
## R2m R2c
## theoretical 0.3232656 0.9128366
## delta 0.3111874 0.8787302Diagnostico do modelo cheio mais plausível
Figura 2.3 Resíduos Quantílicos do modelo cheio plausível
Figura 2.4 Quantile-quantile plot random effects.
Figura 2.5 Predito e observado pelo modelo. Em vermelho linha y=x; em azul um modelo linear por sítio de amostragem.
Figura 2.6 Predito e observado para modelo cheio
Subset=MNEE
Estudo GAMM: smoother type
- Antes de comparar os modelos cheios preciso avaliar smoother type mais para utilizar na GAMM;
- Pelo observado nos dados (figuras 2.1 e 2.2) vou comparar dois smoother type: “tp” e “cr”
- Para isso vou ajuster um smoother para d por Sítio de amostragem com um mesmo parâmetro de penalização, dessa maneira eu considero a variação por sítio de amostragem na variável d sem um elevado custo computacional.
# dados
df_md <- df_resultados %>% filter(MN=="EE") %>% distinct()
df_md <- rbind(mutate(df_md,model_class = "cr"),
mutate(df_md,model_class = "tp"))
df_md$model_class <- factor(df_md$model_class)
# função
f_md <- function(dados){
model_class <- unique(dados$model_class)
if(model_class == "cr"){
md_ <- gam(cbind(n_nRef,100-n_nRef) ~
s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
s(d.z,SiteCode,bs="fs",xt=list(bs="cr"), m=1),
data = dados, family = "binomial")
}else{
md_ <- gam(cbind(n_nRef,100-n_nRef) ~
s(p.z,bs="tp") + s(d.z,bs="tp") + ti(p.z,d.z,bs=c("tp","tp")) +
s(d.z,SiteCode,bs="fs",xt=list(bs="tp"), m=1),
data = dados, family = "binomial")
}
return(md_)
}
registerDoMC(2)
l_nRefEE_estudoGAMM <- dlply(df_md,"model_class",f_md,.parallel = TRUE)Comparação Modelos Cheios
# dados
df_md <- df_resultados %>% filter(MN=="EE") %>%
distinct() %>%
mutate(rep = 1)
df_md <- left_join(x=data.frame(rep=1,
md_class=c("glmm d+d^2|Site",
"glmm d|Site",
"glmm 1|Site",
"gamm d|Site for each",
"gamm d|Site common",
"gamm 1|Site")),
y=df_md,
by="rep")
# função de ajuste
f_md <- function(dados){
md_class <- unique(dados$md_class)
if(md_class == "glmm d+d^2|Site"){
md_ <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (d.z + I(d.z^2)) +
(d.z + I(d.z^2)|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e6)),
na.action = "na.fail")
}else if(md_class == "glmm d|Site"){
md_ <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (d.z + I(d.z^2)) +
(d.z|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e6)),
na.action = "na.fail")
}else if(md_class == "glmm 1|Site"){
md_ <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (d.z + I(d.z^2)) +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e6)),
na.action = "na.fail")
}else if(md_class == "gamm d|Site for each"){
md_ <- gam(cbind(n_nRef,100-n_nRef) ~
s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
s(SiteCode,bs="re") + s(d.z,by=SiteCode,bs="cr",m=1),
data = dados, family = "binomial")
}else if(md_class == "gamm d|Site common"){
md_ <- gam(cbind(n_nRef,100-n_nRef) ~
s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
s(d.z,SiteCode,bs="fs",xt=list(bs="cr"), m=1),
data = dados, family = "binomial")
}else{
md_ <- gam(cbind(n_nRef,100-n_nRef) ~
s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
s(SiteCode,bs="re"),
data = dados, family = "binomial")
}
return(md_)
}registerDoMC(3)
l_nRefEE_mdCheio <- dlply(df_md,"md_class",f_md,.parallel = TRUE)## dAICc df weight
## gamm d|Site for each 0.0 668.986466713628 1
## gamm d|Site common 816.9 851.710244484507 <0.001
## glmm d+d^2|Site 5223.6 15 <0.001
## glmm d|Site 8188.0 12 <0.001
## gamm 1|Site 25233.8 126.840686107377 <0.001
## glmm 1|Site 27131.4 10 <0.001Diagnostico modelo cheio
Figura 2.2.1 gam.check(md_nRefEE1)
##
## Method: UBRE Optimizer: outer newton
## full convergence after 33 iterations.
## Gradient range [-5.585756e-06,5.869852e-06]
## (score 0.6972787 & scale 1).
## eigenvalue range [-1.541502e-11,0.0008330883].
## Model rank = 1063 / 1065
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(p.z) 9.00e+00 1.00e+00 1.34 1.00
## s(d.z) 9.00e+00 8.92e+00 1.04 0.98
## ti(p.z,d.z) 1.60e+01 2.22e+00 0.99 0.56
## s(SiteCode) 1.03e+02 1.01e+02 NA NA
## s(d.z):SiteCodeBAjuss 9.00e+00 3.52e+00 1.04 0.99
## s(d.z):SiteCodeBAlenc4 9.00e+00 6.21e+00 1.04 0.97
## s(d.z):SiteCodeBAuruc 9.00e+00 3.92e+00 1.04 0.98
## s(d.z):SiteCodeESsoor 9.00e+00 7.09e+00 1.04 0.97
## s(d.z):SiteCodeMGipia1 9.00e+00 5.60e+00 1.04 0.98
## s(d.z):SiteCodeMGlavr2 9.00e+00 7.39e+00 1.04 0.96
## s(d.z):SiteCodeMGlavr3 9.00e+00 1.00e+00 1.04 0.98
## s(d.z):SiteCodeMGlavr6 9.00e+00 1.77e+00 1.04 0.98
## s(d.z):SiteCodeMGubera 9.00e+00 8.90e+00 1.04 0.97
## s(d.z):SiteCodeMGuberl1 9.00e+00 3.39e+00 1.04 0.98
## s(d.z):SiteCodeMGuberl3 9.00e+00 2.48e+00 1.04 0.97
## s(d.z):SiteCodeMGuberl5 9.00e+00 8.91e+00 1.04 0.96
## s(d.z):SiteCodeMGuberl6 9.00e+00 5.25e+00 1.04 0.98
## s(d.z):SiteCodeMGuberl7 9.00e+00 4.38e+00 1.04 0.98
## s(d.z):SiteCodeMGvico1 9.00e+00 8.53e+00 1.04 0.97
## s(d.z):SiteCodeMGvico16 9.00e+00 8.15e+00 1.04 0.98
## s(d.z):SiteCodeMGvico3 9.00e+00 7.35e+00 1.04 0.98
## s(d.z):SiteCodeMGvico4 9.00e+00 8.13e+00 1.04 0.96
## s(d.z):SiteCodeMSdour 9.00e+00 1.63e-06 1.04 0.97
## s(d.z):SiteCodePEalia 9.00e+00 6.03e+00 1.04 0.98
## s(d.z):SiteCodePEbrejo 9.00e+00 1.29e+00 1.04 0.99
## s(d.z):SiteCodePEcabo2 9.00e+00 4.49e+00 1.04 0.96
## s(d.z):SiteCodePEcaru1 9.00e+00 5.68e+00 1.04 0.97
## s(d.z):SiteCodePEmata2 9.00e+00 7.83e+00 1.04 0.95
## s(d.z):SiteCodePEsvfer 9.00e+00 8.84e+00 1.04 0.98
## s(d.z):SiteCodePRanto10 9.00e+00 1.00e+00 1.04 0.97
## s(d.z):SiteCodePRanto11 9.00e+00 7.61e+00 1.04 0.97
## s(d.z):SiteCodePRanto12 9.00e+00 3.49e+00 1.04 0.97
## s(d.z):SiteCodePRanto13 9.00e+00 4.82e+00 1.04 0.95
## s(d.z):SiteCodePRanto14 9.00e+00 7.74e+00 1.04 0.99
## s(d.z):SiteCodePRanto15 9.00e+00 3.87e+00 1.04 0.99
## s(d.z):SiteCodePRanto4 9.00e+00 1.00e+00 1.04 0.96
## s(d.z):SiteCodePRanto5 9.00e+00 4.54e+00 1.04 0.98
## s(d.z):SiteCodePRanto6 9.00e+00 3.78e+00 1.04 0.98
## s(d.z):SiteCodePRanto7 9.00e+00 2.43e+00 1.04 0.98
## s(d.z):SiteCodePRanto8 9.00e+00 1.18e+00 1.04 0.98
## s(d.z):SiteCodePRanto9 9.00e+00 1.00e+00 1.04 0.97
## s(d.z):SiteCodePRdiam2 9.00e+00 2.93e+00 1.04 0.97
## s(d.z):SiteCodePRdiam3 9.00e+00 3.51e+00 1.04 0.97
## s(d.z):SiteCodePRdiam4 9.00e+00 6.33e+00 1.04 0.98
## s(d.z):SiteCodePRrico1 9.00e+00 6.55e+00 1.04 0.99
## s(d.z):SiteCodePRrico3 9.00e+00 7.74e+00 1.04 0.98
## s(d.z):SiteCodePRsapo 9.00e+00 8.78e+00 1.04 0.98
## s(d.z):SiteCodePRtele 9.00e+00 7.45e+00 1.04 0.98
## s(d.z):SiteCodePRtiba1 9.00e+00 7.90e+00 1.04 0.97
## s(d.z):SiteCodePRvent 9.00e+00 7.32e+00 1.04 0.98
## s(d.z):SiteCodeRJpnrj 9.00e+00 7.68e+00 1.04 0.95
## s(d.z):SiteCodeRScach1 9.00e+00 6.16e+00 1.04 0.96
## s(d.z):SiteCodeRScach2 9.00e+00 8.10e+00 1.04 0.98
## s(d.z):SiteCodeRScach3 9.00e+00 8.93e+00 1.04 0.97
## s(d.z):SiteCodeRScach4 9.00e+00 4.41e+00 1.04 0.96
## s(d.z):SiteCodeRScama 9.00e+00 7.67e+00 1.04 0.96
## s(d.z):SiteCodeRScamb1 9.00e+00 1.00e+00 1.04 0.97
## s(d.z):SiteCodeRScris 9.00e+00 6.99e+00 1.04 0.97
## s(d.z):SiteCodeRSfaxi 9.00e+00 8.31e+00 1.04 0.98
## s(d.z):SiteCodeRSgaur 9.00e+00 7.26e+00 1.04 0.98
## s(d.z):SiteCodeRSmaqu4 9.00e+00 4.45e+00 1.04 0.98
## s(d.z):SiteCodeRSmorr1 9.00e+00 3.74e+00 1.04 0.98
## s(d.z):SiteCodeRSpet1 9.00e+00 5.52e+00 1.04 0.96
## s(d.z):SiteCodeRSrioz1 9.00e+00 6.84e+00 1.04 0.96
## s(d.z):SiteCodeRSsetape 9.00e+00 7.58e+00 1.04 0.98
## s(d.z):SiteCodeRSsini 9.00e+00 2.56e+00 1.04 0.97
## s(d.z):SiteCodeRSsmar2 9.00e+00 3.14e+00 1.04 0.97
## s(d.z):SiteCodeRSvsol 9.00e+00 7.61e+00 1.04 0.98
## s(d.z):SiteCodeSCarar 9.00e+00 6.02e+00 1.04 0.97
## s(d.z):SiteCodeSCcric1 9.00e+00 7.91e+00 1.04 0.97
## s(d.z):SiteCodeSCilho 9.00e+00 6.28e+00 1.04 0.96
## s(d.z):SiteCodeSCilho2 9.00e+00 4.90e+00 1.04 0.98
## s(d.z):SiteCodeSCorle 9.00e+00 5.29e+00 1.04 0.97
## s(d.z):SiteCodeSCserra4 9.00e+00 1.00e+00 1.04 0.97
## s(d.z):SiteCodeSCside 9.00e+00 1.00e+00 1.04 0.96
## s(d.z):SiteCodeSCside1 9.00e+00 4.10e+00 1.04 0.98
## s(d.z):SiteCodeSCside4 9.00e+00 1.00e+00 1.04 0.99
## s(d.z):SiteCodeSCsjo1 9.00e+00 5.63e+00 1.04 0.98
## s(d.z):SiteCodeSCtimbe1 9.00e+00 2.26e+00 1.04 0.99
## s(d.z):SiteCodeSCvolta1 9.00e+00 7.79e+00 1.04 0.99
## s(d.z):SiteCodeSCvolta2 9.00e+00 4.23e+00 1.04 0.99
## s(d.z):SiteCodeSPcana1 9.00e+00 8.81e+00 1.04 0.96
## s(d.z):SiteCodeSPcara1 9.00e+00 5.15e+00 1.04 0.99
## s(d.z):SiteCodeSPcara3 9.00e+00 8.76e+00 1.04 0.99
## s(d.z):SiteCodeSPeec1 9.00e+00 7.99e+00 1.04 0.99
## s(d.z):SiteCodeSPeec4 9.00e+00 8.34e+00 1.04 0.96
## s(d.z):SiteCodeSPeec5 9.00e+00 7.12e+00 1.04 0.96
## s(d.z):SiteCodeSPigua1 9.00e+00 5.33e+00 1.04 0.98
## s(d.z):SiteCodeSPjabo1 9.00e+00 8.24e+00 1.04 0.97
## s(d.z):SiteCodeSPpecb1 9.00e+00 7.07e+00 1.04 0.97
## s(d.z):SiteCodeSPpei1 9.00e+00 6.97e+00 1.04 0.98
## s(d.z):SiteCodeSPpesm6 9.00e+00 8.54e+00 1.04 0.96
## s(d.z):SiteCodeSPpesmA 9.00e+00 7.99e+00 1.04 0.96
## s(d.z):SiteCodeSPpesmB 9.00e+00 8.05e+00 1.04 0.98
## s(d.z):SiteCodeSPpesmC 9.00e+00 3.27e+00 1.04 0.99
## s(d.z):SiteCodeSPpesmD 9.00e+00 1.00e+00 1.04 0.97
## s(d.z):SiteCodeSPpesmE 9.00e+00 3.89e+00 1.04 0.97
## s(d.z):SiteCodeSPpesmF 9.00e+00 6.67e+00 1.04 0.99
## s(d.z):SiteCodeSPpesmG 9.00e+00 7.21e+00 1.04 0.99
## s(d.z):SiteCodeSPpesmH 9.00e+00 1.83e+00 1.04 0.99
## s(d.z):SiteCodeSPpesmI 9.00e+00 5.21e+00 1.04 0.98
## s(d.z):SiteCodeSPpesmJ 9.00e+00 4.51e+00 1.04 0.95
## s(d.z):SiteCodeSPpesmK 9.00e+00 4.79e+00 1.04 0.96
## s(d.z):SiteCodeSPpesmL 9.00e+00 7.75e+00 1.04 0.99
## s(d.z):SiteCodeSPpesmN 9.00e+00 1.00e+00 1.04 0.98
## s(d.z):SiteCodeSPpesmS 9.00e+00 6.23e+00 1.04 0.99
## s(d.z):SiteCodeSPpesmU 9.00e+00 1.00e+00 1.04 0.99tabela 2.2.1 Summary
##
## Family: binomial
## Link function: logit
##
## Formula:
## cbind(n_nRef, 100 - n_nRef) ~ s(p.z, bs = "cr") + s(d.z, bs = "cr") +
## ti(p.z, d.z, bs = c("cr", "cr")) + s(SiteCode, bs = "re") +
## s(d.z, by = SiteCode, bs = "cr", m = 1)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.985 334.817 0.009 0.993
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(p.z) 1.000e+00 1.000e+00 0.000 0.992501
## s(d.z) 8.917e+00 8.981e+00 184.902 < 2e-16 ***
## ti(p.z,d.z) 2.217e+00 2.550e+00 5.655 0.090208 .
## s(SiteCode) 1.007e+02 1.010e+02 16720.810 < 2e-16 ***
## s(d.z):SiteCodeBAjuss 3.518e+00 4.140e+00 55.864 3.48e-11 ***
## s(d.z):SiteCodeBAlenc4 6.213e+00 6.871e+00 22.646 0.001850 **
## s(d.z):SiteCodeBAuruc 3.918e+00 4.725e+00 16.256 0.004800 **
## s(d.z):SiteCodeESsoor 7.086e+00 7.898e+00 71.275 3.94e-12 ***
## s(d.z):SiteCodeMGipia1 5.600e+00 6.418e+00 19.195 0.005755 **
## s(d.z):SiteCodeMGlavr2 7.393e+00 7.876e+00 101.386 < 2e-16 ***
## s(d.z):SiteCodeMGlavr3 1.000e+00 1.000e+00 0.000 0.999850
## s(d.z):SiteCodeMGlavr6 1.765e+00 1.948e+00 0.122 0.942107
## s(d.z):SiteCodeMGubera 8.901e+00 8.992e+00 374.208 < 2e-16 ***
## s(d.z):SiteCodeMGuberl1 3.389e+00 3.908e+00 23.179 0.000105 ***
## s(d.z):SiteCodeMGuberl3 2.480e+00 3.080e+00 9.264 0.025249 *
## s(d.z):SiteCodeMGuberl5 8.906e+00 8.992e+00 128.865 < 2e-16 ***
## s(d.z):SiteCodeMGuberl6 5.254e+00 5.888e+00 9.589 0.115697
## s(d.z):SiteCodeMGuberl7 4.385e+00 4.918e+00 70.432 1.19e-13 ***
## s(d.z):SiteCodeMGvico1 8.525e+00 8.916e+00 213.180 < 2e-16 ***
## s(d.z):SiteCodeMGvico16 8.150e+00 8.782e+00 118.122 < 2e-16 ***
## s(d.z):SiteCodeMGvico3 7.348e+00 8.071e+00 28.360 0.000512 ***
## s(d.z):SiteCodeMGvico4 8.133e+00 8.698e+00 170.410 < 2e-16 ***
## s(d.z):SiteCodeMSdour 1.628e-06 3.254e-06 0.000 1.000000
## s(d.z):SiteCodePEalia 6.033e+00 6.761e+00 19.821 0.005221 **
## s(d.z):SiteCodePEbrejo 1.294e+00 1.533e+00 0.263 0.782471
## s(d.z):SiteCodePEcabo2 4.492e+00 5.453e+00 10.946 0.064084 .
## s(d.z):SiteCodePEcaru1 5.676e+00 6.434e+00 8.413 0.236611
## s(d.z):SiteCodePEmata2 7.827e+00 8.582e+00 114.401 < 2e-16 ***
## s(d.z):SiteCodePEsvfer 8.837e+00 8.974e+00 1.257 0.998538
## s(d.z):SiteCodePRanto10 1.000e+00 1.001e+00 0.121 0.728520
## s(d.z):SiteCodePRanto11 7.611e+00 8.318e+00 32.710 8.71e-05 ***
## s(d.z):SiteCodePRanto12 3.490e+00 4.185e+00 6.615 0.170366
## s(d.z):SiteCodePRanto13 4.820e+00 5.781e+00 21.908 0.001183 **
## s(d.z):SiteCodePRanto14 7.735e+00 8.369e+00 36.584 1.58e-05 ***
## s(d.z):SiteCodePRanto15 3.869e+00 4.386e+00 6.884 0.215468
## s(d.z):SiteCodePRanto4 1.000e+00 1.000e+00 0.118 0.731777
## s(d.z):SiteCodePRanto5 4.537e+00 5.491e+00 15.817 0.010771 *
## s(d.z):SiteCodePRanto6 3.780e+00 4.453e+00 30.756 7.98e-06 ***
## s(d.z):SiteCodePRanto7 2.431e+00 2.915e+00 13.806 0.006562 **
## s(d.z):SiteCodePRanto8 1.181e+00 1.789e+00 1.340 0.496320
## s(d.z):SiteCodePRanto9 1.000e+00 1.000e+00 0.130 0.718262
## s(d.z):SiteCodePRdiam2 2.929e+00 3.334e+00 7.946 0.070747 .
## s(d.z):SiteCodePRdiam3 3.509e+00 4.042e+00 37.371 1.60e-07 ***
## s(d.z):SiteCodePRdiam4 6.334e+00 7.126e+00 8.110 0.334557
## s(d.z):SiteCodePRrico1 6.553e+00 7.391e+00 98.485 < 2e-16 ***
## s(d.z):SiteCodePRrico3 7.743e+00 8.230e+00 81.085 2.30e-13 ***
## s(d.z):SiteCodePRsapo 8.782e+00 8.967e+00 94.450 < 2e-16 ***
## s(d.z):SiteCodePRtele 7.446e+00 8.185e+00 127.443 < 2e-16 ***
## s(d.z):SiteCodePRtiba1 7.895e+00 8.546e+00 130.425 < 2e-16 ***
## s(d.z):SiteCodePRvent 7.316e+00 7.892e+00 137.352 < 2e-16 ***
## s(d.z):SiteCodeRJpnrj 7.676e+00 8.333e+00 71.624 4.39e-12 ***
## s(d.z):SiteCodeRScach1 6.162e+00 6.781e+00 52.604 7.60e-09 ***
## s(d.z):SiteCodeRScach2 8.105e+00 8.620e+00 156.415 < 2e-16 ***
## s(d.z):SiteCodeRScach3 8.930e+00 8.997e+00 214.992 < 2e-16 ***
## s(d.z):SiteCodeRScach4 4.410e+00 4.998e+00 9.066 0.113431
## s(d.z):SiteCodeRScama 7.671e+00 8.277e+00 17.616 0.028409 *
## s(d.z):SiteCodeRScamb1 1.000e+00 1.000e+00 0.000 0.995301
## s(d.z):SiteCodeRScris 6.991e+00 7.624e+00 45.902 6.73e-07 ***
## s(d.z):SiteCodeRSfaxi 8.312e+00 8.760e+00 172.834 < 2e-16 ***
## s(d.z):SiteCodeRSgaur 7.264e+00 7.919e+00 31.772 9.40e-05 ***
## s(d.z):SiteCodeRSmaqu4 4.452e+00 5.426e+00 44.325 8.46e-08 ***
## s(d.z):SiteCodeRSmorr1 3.741e+00 4.418e+00 7.862 0.111616
## s(d.z):SiteCodeRSpet1 5.521e+00 6.244e+00 50.795 6.03e-09 ***
## s(d.z):SiteCodeRSrioz1 6.841e+00 7.188e+00 0.515 0.999512
## s(d.z):SiteCodeRSsetape 7.577e+00 7.892e+00 14.478 0.074150 .
## s(d.z):SiteCodeRSsini 2.559e+00 3.215e+00 5.872 0.091947 .
## s(d.z):SiteCodeRSsmar2 3.143e+00 3.599e+00 12.399 0.019760 *
## s(d.z):SiteCodeRSvsol 7.609e+00 8.296e+00 43.975 7.49e-07 ***
## s(d.z):SiteCodeSCarar 6.021e+00 6.762e+00 7.108 0.390389
## s(d.z):SiteCodeSCcric1 7.914e+00 8.496e+00 116.890 < 2e-16 ***
## s(d.z):SiteCodeSCilho 6.284e+00 7.056e+00 49.581 1.92e-08 ***
## s(d.z):SiteCodeSCilho2 4.900e+00 5.863e+00 17.084 0.008020 **
## s(d.z):SiteCodeSCorle 5.294e+00 5.952e+00 77.425 1.17e-14 ***
## s(d.z):SiteCodeSCserra4 1.000e+00 1.001e+00 0.000 0.993642
## s(d.z):SiteCodeSCside 1.000e+00 1.001e+00 0.000 0.990357
## s(d.z):SiteCodeSCside1 4.103e+00 4.842e+00 38.398 3.26e-07 ***
## s(d.z):SiteCodeSCside4 1.000e+00 1.000e+00 0.000 0.996727
## s(d.z):SiteCodeSCsjo1 5.627e+00 6.282e+00 27.456 0.000167 ***
## s(d.z):SiteCodeSCtimbe1 2.263e+00 2.833e+00 5.977 0.177634
## s(d.z):SiteCodeSCvolta1 7.792e+00 8.490e+00 102.101 < 2e-16 ***
## s(d.z):SiteCodeSCvolta2 4.227e+00 5.001e+00 18.193 0.002718 **
## s(d.z):SiteCodeSPcana1 8.810e+00 8.981e+00 108.511 < 2e-16 ***
## s(d.z):SiteCodeSPcara1 5.154e+00 6.079e+00 17.558 0.007704 **
## s(d.z):SiteCodeSPcara3 8.763e+00 8.966e+00 126.650 < 2e-16 ***
## s(d.z):SiteCodeSPeec1 7.991e+00 8.609e+00 132.468 < 2e-16 ***
## s(d.z):SiteCodeSPeec4 8.341e+00 8.786e+00 25.417 0.001971 **
## s(d.z):SiteCodeSPeec5 7.123e+00 7.828e+00 45.583 3.12e-07 ***
## s(d.z):SiteCodeSPigua1 5.332e+00 6.161e+00 14.986 0.025390 *
## s(d.z):SiteCodeSPjabo1 8.242e+00 8.718e+00 255.351 < 2e-16 ***
## s(d.z):SiteCodeSPpecb1 7.066e+00 7.735e+00 50.802 2.18e-08 ***
## s(d.z):SiteCodeSPpei1 6.966e+00 7.873e+00 74.614 6.27e-13 ***
## s(d.z):SiteCodeSPpesm6 8.537e+00 8.898e+00 114.199 < 2e-16 ***
## s(d.z):SiteCodeSPpesmA 7.988e+00 8.645e+00 41.836 3.61e-06 ***
## s(d.z):SiteCodeSPpesmB 8.046e+00 8.598e+00 14.627 0.139066
## s(d.z):SiteCodeSPpesmC 3.271e+00 3.928e+00 4.808 0.239948
## s(d.z):SiteCodeSPpesmD 1.000e+00 1.001e+00 0.000 0.993652
## s(d.z):SiteCodeSPpesmE 3.895e+00 4.591e+00 6.616 0.182925
## s(d.z):SiteCodeSPpesmF 6.674e+00 7.643e+00 11.998 0.128377
## s(d.z):SiteCodeSPpesmG 7.208e+00 7.967e+00 20.811 0.009023 **
## s(d.z):SiteCodeSPpesmH 1.825e+00 2.285e+00 1.549 0.460585
## s(d.z):SiteCodeSPpesmI 5.206e+00 6.062e+00 20.324 0.002539 **
## s(d.z):SiteCodeSPpesmJ 4.505e+00 5.275e+00 8.113 0.173980
## s(d.z):SiteCodeSPpesmK 4.786e+00 5.615e+00 27.558 6.19e-05 ***
## s(d.z):SiteCodeSPpesmL 7.750e+00 8.423e+00 38.378 9.06e-06 ***
## s(d.z):SiteCodeSPpesmN 1.000e+00 1.001e+00 0.007 0.932461
## s(d.z):SiteCodeSPpesmS 6.230e+00 7.199e+00 8.561 0.292272
## s(d.z):SiteCodeSPpesmU 1.000e+00 1.000e+00 0.000 0.985334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 1063/1065
## R-sq.(adj) = 0.986 Deviance explained = 98.3%
## UBRE = 0.69728 Scale est. = 1 n = 2060Figura 2.2.2 Observado e predito pelo modelo cheio mais plausível
Figura 2.2.3 Diagnostico: avaliação do ajuste pelas preditoras lineares
v_sitesOut_nRefEE <- filter(df_plot,.fitted < -300 | .fitted > 380 | .cooksd > 2000) %>%
.$SiteCode %>% unique
Comparação Modelo Cheio sem outliers
# sites outlier
v_sitesOut_nRefEE <- filter(df_plot,.fitted < -300 | .fitted > 380 | .cooksd > 2000) %>%
.$SiteCode %>% unique
# dados
df_md <- df_resultados %>% filter(MN=="EE" & !(SiteCode %in% v_sitesOut_nRefEE)) %>%
distinct() %>%
mutate(rep = 1)
df_md <- left_join(x=data.frame(rep=1,
md_class=c("glmm d+d^2|Site",
"glmm d|Site",
"glmm 1|Site",
"gamm d|Site for each",
"gamm d|Site common",
"gamm 1|Site")),
y=df_md,
by="rep")
# ajuste modelos
registerDoMC(3)
l_nRefEE_mdCheio__sOut <- dlply(df_md,"md_class",f_md,.parallel = TRUE)## dAICc df weight
## gamm d|Site for each 0.0 644.383383100452 1
## gamm d|Site common 778.8 826.614698460209 <0.001
## glmm d+d^2|Site 5111.4 15 <0.001
## glmm d|Site 8063.3 12 <0.001
## gamm 1|Site 23983.9 123.849903818296 <0.001
## glmm 1|Site 26005.3 10 <0.001Diagnosticos modelo cheio sem outliers
Figura 2.2.4 Appraise(md_nRefEE_sOut)
Figura 2.2.5 Diagnostico: avaliação do ajuste pelas preditoras lineares
Subset=MNEI
Estudo GAMM: smoother type
- Antes de comparar os modelos cheios preciso avaliar smoother type mais para utilizar na GAMM;
- Pelo observado nos dados (figuras 2.1 e 2.2) vou comparar dois smoother type: “tp” e “cr”
- Para isso vou ajuster um smoother para d por Sítio de amostragem com um mesmo parâmetro de penalização, dessa maneira eu considero a variação por sítio de amostragem na variável d sem um elevado custo computacional.
# dados
df_md <- df_resultados %>% filter(MN=="EI") %>% distinct()
df_md <- rbind(mutate(df_md,model_class = "cr"),
mutate(df_md,model_class = "tp"))
df_md$model_class <- factor(df_md$model_class)
# função
f_md <- function(dados){
model_class <- unique(dados$model_class)
if(model_class == "cr"){
md_ <- gam(cbind(n_nRef,100-n_nRef) ~
s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
s(d.z,SiteCode,bs="fs",xt=list(bs="cr"), m=1),
data = dados, family = "binomial")
}else{
md_ <- gam(cbind(n_nRef,100-n_nRef) ~
s(p.z,bs="tp") + s(d.z,bs="tp") + ti(p.z,d.z,bs=c("tp","tp")) +
s(d.z,SiteCode,bs="fs",xt=list(bs="tp"), m=1),
data = dados, family = "binomial")
}
return(md_)
}
registerDoMC(2)
l_nRefEI_estudoGAMM <- dlply(df_md,"model_class",f_md,.parallel = TRUE)## dAICc df weight
## cr 0.0 647.199595731505 1
## tp 443.6 758.152972168071 <0.001Comparação Modelos Cheios
# dados
df_md <- df_resultados %>% filter(MN=="EI") %>%
distinct() %>%
mutate(rep = 1)
df_md <- left_join(x=data.frame(rep=1,
md_class=c("glmm d+d^2|Site",
"glmm d|Site",
"glmm 1|Site",
"gamm d|Site for each",
"gamm d|Site common",
"gamm 1|Site")),
y=df_md,
by="rep")
# função de ajuste
registerDoMC(3)
l_nRefEI_mdCheio <- dlply(df_md,"md_class",f_md,.parallel = TRUE)## dAICc df weight
## gamm d|Site for each 0.0 556.021191967948 1
## gamm d|Site common 374.1 647.199595731505 <0.001
## glmm d+d^2|Site 3256.4 15 <0.001
## glmm d|Site 11864.4 12 <0.001
## gamm 1|Site 57181.1 125.994597023309 <0.001
## glmm 1|Site 63995.3 10 <0.001Diagnostico modelo cheio
Figura 2.3.1 appraise(md_nRefEI)
##
## Method: UBRE Optimizer: outer newton
## full convergence after 50 iterations.
## Gradient range [-2.37292e-06,8.788319e-07]
## (score 0.1963635 & scale 1).
## eigenvalue range [-7.886026e-07,0.001085589].
## Model rank = 1063 / 1065
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(p.z) 9.00e+00 1.00e+00 1.25 1.000
## s(d.z) 9.00e+00 7.25e+00 0.96 0.020 *
## ti(p.z,d.z) 1.60e+01 1.29e+01 1.00 0.640
## s(SiteCode) 1.03e+02 1.01e+02 NA NA
## s(d.z):SiteCodeBAjuss 9.00e+00 2.88e+00 0.96 0.040 *
## s(d.z):SiteCodeBAlenc4 9.00e+00 3.45e+00 0.96 0.025 *
## s(d.z):SiteCodeBAuruc 9.00e+00 1.00e+00 0.96 0.045 *
## s(d.z):SiteCodeESsoor 9.00e+00 2.52e+00 0.96 0.015 *
## s(d.z):SiteCodeMGipia1 9.00e+00 3.53e+00 0.96 0.065 .
## s(d.z):SiteCodeMGlavr2 9.00e+00 8.56e+00 0.96 0.035 *
## s(d.z):SiteCodeMGlavr3 9.00e+00 5.20e-05 0.96 0.020 *
## s(d.z):SiteCodeMGlavr6 9.00e+00 7.22e+00 0.96 0.015 *
## s(d.z):SiteCodeMGubera 9.00e+00 6.97e+00 0.96 0.040 *
## s(d.z):SiteCodeMGuberl1 9.00e+00 4.87e+00 0.96 0.050 *
## s(d.z):SiteCodeMGuberl3 9.00e+00 4.06e+00 0.96 0.015 *
## s(d.z):SiteCodeMGuberl5 9.00e+00 5.98e+00 0.96 0.020 *
## s(d.z):SiteCodeMGuberl6 9.00e+00 4.05e+00 0.96 0.025 *
## s(d.z):SiteCodeMGuberl7 9.00e+00 3.50e+00 0.96 0.045 *
## s(d.z):SiteCodeMGvico1 9.00e+00 1.00e+00 0.96 0.025 *
## s(d.z):SiteCodeMGvico16 9.00e+00 1.00e+00 0.96 0.045 *
## s(d.z):SiteCodeMGvico3 9.00e+00 1.07e+00 0.96 0.050 *
## s(d.z):SiteCodeMGvico4 9.00e+00 1.00e+00 0.96 0.025 *
## s(d.z):SiteCodeMSdour 9.00e+00 2.53e+00 0.96 0.035 *
## s(d.z):SiteCodePEalia 9.00e+00 2.41e+00 0.96 0.045 *
## s(d.z):SiteCodePEbrejo 9.00e+00 5.62e+00 0.96 0.015 *
## s(d.z):SiteCodePEcabo2 9.00e+00 4.30e+00 0.96 0.045 *
## s(d.z):SiteCodePEcaru1 9.00e+00 5.14e+00 0.96 0.070 .
## s(d.z):SiteCodePEmata2 9.00e+00 1.00e+00 0.96 0.045 *
## s(d.z):SiteCodePEsvfer 9.00e+00 8.26e+00 0.96 0.050 *
## s(d.z):SiteCodePRanto10 9.00e+00 6.33e+00 0.96 0.045 *
## s(d.z):SiteCodePRanto11 9.00e+00 3.68e+00 0.96 0.055 .
## s(d.z):SiteCodePRanto12 9.00e+00 7.30e+00 0.96 0.040 *
## s(d.z):SiteCodePRanto13 9.00e+00 1.00e+00 0.96 0.025 *
## s(d.z):SiteCodePRanto14 9.00e+00 5.90e+00 0.96 0.040 *
## s(d.z):SiteCodePRanto15 9.00e+00 2.85e+00 0.96 0.040 *
## s(d.z):SiteCodePRanto4 9.00e+00 5.53e+00 0.96 0.040 *
## s(d.z):SiteCodePRanto5 9.00e+00 1.00e+00 0.96 0.030 *
## s(d.z):SiteCodePRanto6 9.00e+00 6.85e+00 0.96 0.045 *
## s(d.z):SiteCodePRanto7 9.00e+00 5.64e+00 0.96 0.050 *
## s(d.z):SiteCodePRanto8 9.00e+00 1.00e+00 0.96 0.015 *
## s(d.z):SiteCodePRanto9 9.00e+00 5.01e+00 0.96 0.035 *
## s(d.z):SiteCodePRdiam2 9.00e+00 1.23e+00 0.96 0.020 *
## s(d.z):SiteCodePRdiam3 9.00e+00 1.02e+00 0.96 0.040 *
## s(d.z):SiteCodePRdiam4 9.00e+00 2.88e+00 0.96 0.020 *
## s(d.z):SiteCodePRrico1 9.00e+00 4.73e+00 0.96 0.025 *
## s(d.z):SiteCodePRrico3 9.00e+00 1.00e+00 0.96 0.045 *
## s(d.z):SiteCodePRsapo 9.00e+00 6.54e+00 0.96 0.050 *
## s(d.z):SiteCodePRtele 9.00e+00 5.02e+00 0.96 0.030 *
## s(d.z):SiteCodePRtiba1 9.00e+00 2.65e+00 0.96 0.045 *
## s(d.z):SiteCodePRvent 9.00e+00 5.25e+00 0.96 0.035 *
## s(d.z):SiteCodeRJpnrj 9.00e+00 1.00e+00 0.96 0.070 .
## s(d.z):SiteCodeRScach1 9.00e+00 3.56e+00 0.96 0.040 *
## s(d.z):SiteCodeRScach2 9.00e+00 2.39e+00 0.96 0.045 *
## s(d.z):SiteCodeRScach3 9.00e+00 2.83e+00 0.96 0.015 *
## s(d.z):SiteCodeRScach4 9.00e+00 1.00e+00 0.96 0.060 .
## s(d.z):SiteCodeRScama 9.00e+00 7.58e+00 0.96 0.015 *
## s(d.z):SiteCodeRScamb1 9.00e+00 3.01e+00 0.96 0.030 *
## s(d.z):SiteCodeRScris 9.00e+00 5.40e+00 0.96 0.065 .
## s(d.z):SiteCodeRSfaxi 9.00e+00 4.59e+00 0.96 0.035 *
## s(d.z):SiteCodeRSgaur 9.00e+00 4.13e+00 0.96 0.020 *
## s(d.z):SiteCodeRSmaqu4 9.00e+00 7.67e+00 0.96 0.040 *
## s(d.z):SiteCodeRSmorr1 9.00e+00 6.91e+00 0.96 0.035 *
## s(d.z):SiteCodeRSpet1 9.00e+00 5.21e+00 0.96 0.055 .
## s(d.z):SiteCodeRSrioz1 9.00e+00 5.08e+00 0.96 0.040 *
## s(d.z):SiteCodeRSsetape 9.00e+00 5.86e+00 0.96 0.030 *
## s(d.z):SiteCodeRSsini 9.00e+00 6.37e+00 0.96 0.035 *
## s(d.z):SiteCodeRSsmar2 9.00e+00 8.17e+00 0.96 0.055 .
## s(d.z):SiteCodeRSvsol 9.00e+00 2.73e+00 0.96 0.010 **
## s(d.z):SiteCodeSCarar 9.00e+00 1.51e+00 0.96 0.030 *
## s(d.z):SiteCodeSCcric1 9.00e+00 3.16e+00 0.96 0.025 *
## s(d.z):SiteCodeSCilho 9.00e+00 2.98e+00 0.96 0.045 *
## s(d.z):SiteCodeSCilho2 9.00e+00 4.10e+00 0.96 0.055 .
## s(d.z):SiteCodeSCorle 9.00e+00 5.77e+00 0.96 0.040 *
## s(d.z):SiteCodeSCserra4 9.00e+00 6.62e+00 0.96 0.045 *
## s(d.z):SiteCodeSCside 9.00e+00 3.28e+00 0.96 0.045 *
## s(d.z):SiteCodeSCside1 9.00e+00 2.71e+00 0.96 0.025 *
## s(d.z):SiteCodeSCside4 9.00e+00 1.00e+00 0.96 0.020 *
## s(d.z):SiteCodeSCsjo1 9.00e+00 4.88e+00 0.96 0.060 .
## s(d.z):SiteCodeSCtimbe1 9.00e+00 6.09e+00 0.96 0.035 *
## s(d.z):SiteCodeSCvolta1 9.00e+00 4.50e+00 0.96 0.030 *
## s(d.z):SiteCodeSCvolta2 9.00e+00 5.90e+00 0.96 0.020 *
## s(d.z):SiteCodeSPcana1 9.00e+00 5.20e+00 0.96 0.030 *
## s(d.z):SiteCodeSPcara1 9.00e+00 7.14e+00 0.96 0.030 *
## s(d.z):SiteCodeSPcara3 9.00e+00 7.60e+00 0.96 0.045 *
## s(d.z):SiteCodeSPeec1 9.00e+00 1.00e+00 0.96 0.055 .
## s(d.z):SiteCodeSPeec4 9.00e+00 6.21e+00 0.96 0.035 *
## s(d.z):SiteCodeSPeec5 9.00e+00 1.59e+00 0.96 0.045 *
## s(d.z):SiteCodeSPigua1 9.00e+00 1.00e+00 0.96 0.055 .
## s(d.z):SiteCodeSPjabo1 9.00e+00 6.59e+00 0.96 0.025 *
## s(d.z):SiteCodeSPpecb1 9.00e+00 1.93e+00 0.96 0.020 *
## s(d.z):SiteCodeSPpei1 9.00e+00 5.09e+00 0.96 0.050 *
## s(d.z):SiteCodeSPpesm6 9.00e+00 3.56e+00 0.96 0.035 *
## s(d.z):SiteCodeSPpesmA 9.00e+00 5.78e+00 0.96 0.025 *
## s(d.z):SiteCodeSPpesmB 9.00e+00 3.44e+00 0.96 0.015 *
## s(d.z):SiteCodeSPpesmC 9.00e+00 3.10e+00 0.96 0.020 *
## s(d.z):SiteCodeSPpesmD 9.00e+00 3.66e+00 0.96 0.030 *
## s(d.z):SiteCodeSPpesmE 9.00e+00 4.04e+00 0.96 0.035 *
## s(d.z):SiteCodeSPpesmF 9.00e+00 3.50e+00 0.96 0.065 .
## s(d.z):SiteCodeSPpesmG 9.00e+00 5.09e+00 0.96 0.045 *
## s(d.z):SiteCodeSPpesmH 9.00e+00 5.30e+00 0.96 0.045 *
## s(d.z):SiteCodeSPpesmI 9.00e+00 6.60e+00 0.96 0.035 *
## s(d.z):SiteCodeSPpesmJ 9.00e+00 6.83e+00 0.96 0.045 *
## s(d.z):SiteCodeSPpesmK 9.00e+00 4.66e+00 0.96 0.050 *
## s(d.z):SiteCodeSPpesmL 9.00e+00 5.60e+00 0.96 0.025 *
## s(d.z):SiteCodeSPpesmN 9.00e+00 4.81e+00 0.96 0.025 *
## s(d.z):SiteCodeSPpesmS 9.00e+00 5.19e+00 0.96 0.025 *
## s(d.z):SiteCodeSPpesmU 9.00e+00 8.32e+00 0.96 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Family: binomial
## Link function: logit
##
## Formula:
## cbind(n_nRef, 100 - n_nRef) ~ s(p.z, bs = "cr") + s(d.z, bs = "cr") +
## ti(p.z, d.z, bs = c("cr", "cr")) + s(SiteCode, bs = "re") +
## s(d.z, by = SiteCode, bs = "cr", m = 1)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5485 10.6261 -0.052 0.959
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(p.z) 1.000e+00 1.000e+00 0.001 0.969635
## s(d.z) 7.249e+00 7.995e+00 599.870 < 2e-16 ***
## ti(p.z,d.z) 1.288e+01 1.296e+01 1918.523 < 2e-16 ***
## s(SiteCode) 1.008e+02 1.010e+02 12589.462 < 2e-16 ***
## s(d.z):SiteCodeBAjuss 2.885e+00 3.083e+00 7.211 0.043988 *
## s(d.z):SiteCodeBAlenc4 3.450e+00 4.178e+00 14.025 0.008454 **
## s(d.z):SiteCodeBAuruc 1.000e+00 1.000e+00 2.169 0.140871
## s(d.z):SiteCodeESsoor 2.520e+00 2.947e+00 11.587 0.016453 *
## s(d.z):SiteCodeMGipia1 3.526e+00 4.180e+00 10.011 0.042309 *
## s(d.z):SiteCodeMGlavr2 8.565e+00 8.914e+00 33.426 0.000297 ***
## s(d.z):SiteCodeMGlavr3 5.202e-05 9.743e-05 0.000 0.999819
## s(d.z):SiteCodeMGlavr6 7.224e+00 8.001e+00 19.314 0.013379 *
## s(d.z):SiteCodeMGubera 6.967e+00 7.543e+00 88.994 1.27e-15 ***
## s(d.z):SiteCodeMGuberl1 4.869e+00 5.678e+00 152.128 < 2e-16 ***
## s(d.z):SiteCodeMGuberl3 4.057e+00 4.755e+00 50.761 1.09e-09 ***
## s(d.z):SiteCodeMGuberl5 5.977e+00 6.744e+00 41.273 7.68e-07 ***
## s(d.z):SiteCodeMGuberl6 4.053e+00 4.746e+00 69.693 1.86e-13 ***
## s(d.z):SiteCodeMGuberl7 3.503e+00 4.176e+00 22.220 0.000249 ***
## s(d.z):SiteCodeMGvico1 1.000e+00 1.000e+00 10.331 0.001310 **
## s(d.z):SiteCodeMGvico16 1.000e+00 1.000e+00 10.693 0.001075 **
## s(d.z):SiteCodeMGvico3 1.066e+00 1.127e+00 11.695 0.001110 **
## s(d.z):SiteCodeMGvico4 1.002e+00 1.004e+00 9.709 0.001845 **
## s(d.z):SiteCodeMSdour 2.529e+00 3.209e+00 16.050 0.001154 **
## s(d.z):SiteCodePEalia 2.409e+00 2.912e+00 24.967 0.000137 ***
## s(d.z):SiteCodePEbrejo 5.618e+00 6.529e+00 102.805 < 2e-16 ***
## s(d.z):SiteCodePEcabo2 4.296e+00 5.148e+00 96.974 < 2e-16 ***
## s(d.z):SiteCodePEcaru1 5.137e+00 5.974e+00 125.181 < 2e-16 ***
## s(d.z):SiteCodePEmata2 1.000e+00 1.000e+00 2.394 0.121837
## s(d.z):SiteCodePEsvfer 8.264e+00 8.806e+00 162.156 < 2e-16 ***
## s(d.z):SiteCodePRanto10 6.330e+00 7.239e+00 126.457 < 2e-16 ***
## s(d.z):SiteCodePRanto11 3.679e+00 4.402e+00 27.440 6.02e-05 ***
## s(d.z):SiteCodePRanto12 7.299e+00 8.174e+00 60.495 4.90e-10 ***
## s(d.z):SiteCodePRanto13 1.001e+00 1.002e+00 1.105 0.294043
## s(d.z):SiteCodePRanto14 5.902e+00 7.076e+00 37.792 3.99e-06 ***
## s(d.z):SiteCodePRanto15 2.853e+00 3.519e+00 5.771 0.315062
## s(d.z):SiteCodePRanto4 5.527e+00 6.364e+00 66.910 3.70e-12 ***
## s(d.z):SiteCodePRanto5 1.000e+00 1.001e+00 1.044 0.307079
## s(d.z):SiteCodePRanto6 6.845e+00 7.812e+00 9.459 0.291539
## s(d.z):SiteCodePRanto7 5.640e+00 6.689e+00 359.758 < 2e-16 ***
## s(d.z):SiteCodePRanto8 1.004e+00 1.009e+00 1.117 0.291522
## s(d.z):SiteCodePRanto9 5.009e+00 5.897e+00 144.199 < 2e-16 ***
## s(d.z):SiteCodePRdiam2 1.226e+00 1.418e+00 10.480 0.005616 **
## s(d.z):SiteCodePRdiam3 1.015e+00 1.030e+00 8.877 0.002987 **
## s(d.z):SiteCodePRdiam4 2.878e+00 3.612e+00 13.345 0.010096 *
## s(d.z):SiteCodePRrico1 4.728e+00 5.406e+00 14.379 0.019964 *
## s(d.z):SiteCodePRrico3 1.000e+00 1.000e+00 7.989 0.004707 **
## s(d.z):SiteCodePRsapo 6.541e+00 7.501e+00 330.201 < 2e-16 ***
## s(d.z):SiteCodePRtele 5.021e+00 5.832e+00 17.533 0.005833 **
## s(d.z):SiteCodePRtiba1 2.646e+00 3.260e+00 10.179 0.019416 *
## s(d.z):SiteCodePRvent 5.248e+00 6.146e+00 94.939 < 2e-16 ***
## s(d.z):SiteCodeRJpnrj 1.000e+00 1.000e+00 7.511 0.006136 **
## s(d.z):SiteCodeRScach1 3.557e+00 4.338e+00 34.071 1.13e-06 ***
## s(d.z):SiteCodeRScach2 2.389e+00 3.055e+00 12.198 0.006591 **
## s(d.z):SiteCodeRScach3 2.833e+00 3.439e+00 63.362 6.16e-13 ***
## s(d.z):SiteCodeRScach4 1.000e+00 1.000e+00 9.436 0.002128 **
## s(d.z):SiteCodeRScama 7.584e+00 8.323e+00 209.383 < 2e-16 ***
## s(d.z):SiteCodeRScamb1 3.008e+00 3.665e+00 52.535 7.17e-11 ***
## s(d.z):SiteCodeRScris 5.399e+00 6.118e+00 58.011 1.41e-10 ***
## s(d.z):SiteCodeRSfaxi 4.590e+00 5.314e+00 307.812 < 2e-16 ***
## s(d.z):SiteCodeRSgaur 4.128e+00 4.885e+00 286.344 < 2e-16 ***
## s(d.z):SiteCodeRSmaqu4 7.666e+00 8.415e+00 45.417 4.66e-07 ***
## s(d.z):SiteCodeRSmorr1 6.906e+00 7.821e+00 293.708 < 2e-16 ***
## s(d.z):SiteCodeRSpet1 5.209e+00 6.231e+00 10.955 0.074749 .
## s(d.z):SiteCodeRSrioz1 5.080e+00 5.964e+00 31.943 1.92e-05 ***
## s(d.z):SiteCodeRSsetape 5.856e+00 6.568e+00 167.875 < 2e-16 ***
## s(d.z):SiteCodeRSsini 6.372e+00 7.146e+00 72.276 6.73e-13 ***
## s(d.z):SiteCodeRSsmar2 8.174e+00 8.716e+00 244.028 < 2e-16 ***
## s(d.z):SiteCodeRSvsol 2.727e+00 3.283e+00 51.506 9.54e-10 ***
## s(d.z):SiteCodeSCarar 1.505e+00 1.871e+00 13.075 0.003705 **
## s(d.z):SiteCodeSCcric1 3.163e+00 3.925e+00 17.572 0.001897 **
## s(d.z):SiteCodeSCilho 2.976e+00 3.672e+00 23.916 0.000124 ***
## s(d.z):SiteCodeSCilho2 4.097e+00 4.893e+00 51.290 6.65e-10 ***
## s(d.z):SiteCodeSCorle 5.769e+00 6.826e+00 211.219 < 2e-16 ***
## s(d.z):SiteCodeSCserra4 6.620e+00 7.540e+00 81.477 1.59e-14 ***
## s(d.z):SiteCodeSCside 3.278e+00 4.059e+00 16.239 0.002868 **
## s(d.z):SiteCodeSCside1 2.710e+00 3.297e+00 76.515 5.85e-16 ***
## s(d.z):SiteCodeSCside4 1.000e+00 1.000e+00 0.072 0.788352
## s(d.z):SiteCodeSCsjo1 4.882e+00 5.747e+00 226.634 < 2e-16 ***
## s(d.z):SiteCodeSCtimbe1 6.090e+00 7.086e+00 301.212 < 2e-16 ***
## s(d.z):SiteCodeSCvolta1 4.497e+00 5.367e+00 144.630 < 2e-16 ***
## s(d.z):SiteCodeSCvolta2 5.898e+00 6.835e+00 164.919 < 2e-16 ***
## s(d.z):SiteCodeSPcana1 5.196e+00 6.341e+00 19.944 0.003461 **
## s(d.z):SiteCodeSPcara1 7.135e+00 8.026e+00 326.287 < 2e-16 ***
## s(d.z):SiteCodeSPcara3 7.604e+00 8.432e+00 312.611 < 2e-16 ***
## s(d.z):SiteCodeSPeec1 1.001e+00 1.002e+00 11.864 0.000575 ***
## s(d.z):SiteCodeSPeec4 6.214e+00 7.056e+00 15.967 0.026046 *
## s(d.z):SiteCodeSPeec5 1.585e+00 1.936e+00 9.642 0.005023 **
## s(d.z):SiteCodeSPigua1 1.000e+00 1.000e+00 3.236 0.072067 .
## s(d.z):SiteCodeSPjabo1 6.590e+00 7.380e+00 20.934 0.006802 **
## s(d.z):SiteCodeSPpecb1 1.930e+00 2.544e+00 19.084 0.002274 **
## s(d.z):SiteCodeSPpei1 5.087e+00 6.116e+00 9.397 0.168178
## s(d.z):SiteCodeSPpesm6 3.563e+00 4.371e+00 155.077 < 2e-16 ***
## s(d.z):SiteCodeSPpesmA 5.780e+00 6.802e+00 129.630 < 2e-16 ***
## s(d.z):SiteCodeSPpesmB 3.443e+00 4.133e+00 89.219 < 2e-16 ***
## s(d.z):SiteCodeSPpesmC 3.098e+00 3.763e+00 26.642 2.77e-05 ***
## s(d.z):SiteCodeSPpesmD 3.661e+00 4.266e+00 117.521 < 2e-16 ***
## s(d.z):SiteCodeSPpesmE 4.039e+00 4.853e+00 39.450 2.11e-07 ***
## s(d.z):SiteCodeSPpesmF 3.499e+00 4.188e+00 36.301 3.88e-07 ***
## s(d.z):SiteCodeSPpesmG 5.087e+00 6.079e+00 14.182 0.027445 *
## s(d.z):SiteCodeSPpesmH 5.305e+00 6.211e+00 16.330 0.013242 *
## s(d.z):SiteCodeSPpesmI 6.596e+00 7.549e+00 64.524 7.11e-11 ***
## s(d.z):SiteCodeSPpesmJ 6.829e+00 7.609e+00 66.533 2.29e-11 ***
## s(d.z):SiteCodeSPpesmK 4.660e+00 5.463e+00 53.814 7.52e-10 ***
## s(d.z):SiteCodeSPpesmL 5.601e+00 6.492e+00 86.200 7.25e-16 ***
## s(d.z):SiteCodeSPpesmN 4.807e+00 5.718e+00 37.147 1.12e-06 ***
## s(d.z):SiteCodeSPpesmS 5.188e+00 6.096e+00 281.815 < 2e-16 ***
## s(d.z):SiteCodeSPpesmU 8.322e+00 8.808e+00 242.182 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 1063/1065
## R-sq.(adj) = 0.992 Deviance explained = 99%
## UBRE = 0.19636 Scale est. = 1 n = 2060Figura 2.3.2 Predito e observado por sítio
Figura 2.3.3 Diagnostico: avaliação do ajuste pelas preditoras lineares
2.4 Conjunto Completo dos Dados
Dados
The following code was adapted from Gavin Simpson (2016) blog post [acessado em abril de 2020]:
https://fromthebottomoftheheap.net/2016/12/15/simultaneous-interval-revisited/
# new data
df_pred0 <- expand.grid(SiteCode=levels(df_resultados$SiteCode)[1],
d.z=seq(min(df_resultados$d.z),max(df_resultados$d.z),length=40),
p.z=seq(min(df_resultados$p.z),max(df_resultados$p.z),length=200))
l_mdNref <- list()
l_mdNref[[1]] <- l_nRefEE_mdCheio[["gamm d|Site for each"]]
l_mdNref[[2]] <- l_nRefEI_mdCheio[["gamm d|Site for each"]]
names(l_mdNref) <- c("EE","EI")
# function
f_predict_PostDist <- function(gamm,df=df_pred0,N=1e4){
#objetos em comum
beta <- coef(gamm)
Vb <- vcov(gamm)
Xp <- predict(gamm,df,type="lpmatrix")
pred <- predict.gam(gamm,df,se.fit = TRUE)
se.fit <- pred$se.fit
pred.response <- predict.gam(gamm,newdata = df,type = "response",se.fit = TRUE)
### estabelece um fator de correção para se.fit
BUdiff <- MASS::mvrnorm(N,mu=rep(0,nrow(Vb)),Sigma=Vb)
simDev <- Xp %*% t(BUdiff)
absDev <- abs(sweep(simDev, 1, se.fit, FUN = "/"))
masd <- apply(absDev, 2L, max)
crit <- quantile(masd,prob=0.95,type=8)
### simula todos os intervalos de confinça
### e então calcula os quantis para cada ponto predito
sims <- MASS::mvrnorm(N,mu=beta,Sigma=Vb)
fits <- Xp %*% t(sims)
q_fits <- t(apply(fits,1,quantile,prob=c(.025,0.5,0.975)))
## output
df_return <- cbind(df,
data.frame(fit = pred$fit,
se.fit = pred$se.fit,
fitResponse = pred.response$fit,
se.fitResponse = pred.response$se.fit,
avg_fit = apply(fits,1,mean),
median_fit = q_fits[,2],
q_.025 = q_fits[,1],
q_.975 = q_fits[,3])
)
df_return %<>%
mutate(upper.fit = fit + (crit*se.fit),
lower.fit = fit - (crit*se.fit)) %>%
select(SiteCode,d.z,p.z,fitResponse,se.fitResponse,fit,se.fit,upper.fit,lower.fit,avg_fit,median_fit,q_.025,q_.975)
return(df_return)
}
df_plot <- ldply(l_mdNref,f_predict_PostDist,.id="MN")Figura 2.4.1 Comparação dos métodos de obter os resultados. Na primeira coluna comparação
#
df_plot %<>% mutate(upper.fit_link = fit + (2*se.fit),
lower.fit_link = fit + (2*se.fit))
df_plot1 <- as.data.frame(apply(df_plot[,7:16],2,function(X) arm::invlogit(X)*100))
names(df_plot1) <- sapply(names(df_plot1), function(x) paste0(x,"_response"))
df_plot %<>% cbind(.,df_plot1)
#
l_p <- list()
l_p[[1]] <- ggplot(df_plot,aes(x=fitResponse,y=fit_response)) +
geom_point(alpha=0.3) + labs(title = "mean:invLink(link) ~ response")
l_p[[2]] <- ggplot(df_plot,aes(x=se.fitResponse,y=se.fit_response)) +
geom_point(alpha=0.3) + labs(title = "se:invLink(link) ~ response")
l_p[[3]] <- ggplot(df_plot,aes(x=fit_response,y=median_fit_response)) +
geom_point(alpha=0.3) + labs(title = "median posterior ~ link")
l_p[[4]] <- ggplot(df_plot,aes(x=upper.fit_response,y=q_.025_response)) +
geom_point(alpha=0.3) + labs(title = " 97,5% posterior ~ fit + crit*se.fit")
l_p[[5]] <- ggplot(df_plot,aes(x=lower.fit_response,y=q_.975_response)) +
geom_point(alpha=0.3) + labs(title = "2,5% posterior ~ fit - crit*se.fit")
# do.call("grid.arrange",c(l_p,ncol=3))
grid.arrange(l_p[[1]],l_p[[2]],l_p[[3]],l_p[[4]],l_p[[5]],
layout_matrix=rbind(c(1,3,NA),
c(2,4,5)),
top="response scale")